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Global atmospheric carbon dioxide (C02) measurements for the NASA Active Sensing 
of C02 Emissions over Nights, Days, and Seasons (ASCENDS) space mission are 
critical for improving our understanding of global CO 2 sources and sinks. Advanced 
Intensity- Modulated Continuous-Wave (IM-CW) lidar techniques are investigated as a 
means of facilitating C02 measurements from space to meet the ASCENDS 
measurement requirements. In recent numerical, laboratory and flight experiments we 
have successfully used the Binary Phase Shift Keying (BPSK) modulation technique to 
uniquely discriminate surface lidar returns from intermediate aerosol and cloud 
contamination. We demonstrate the utility of BPSK to eliminate sidelobes in the range 
profile as a means of making Integrated Path Differential Absorption (IPDA) column 
C02 measurements in the presence of optically thin clouds, thereby eliminating the need 
to correct for sidelobe bias errors caused by the clouds. Furthermore, high accuracy and 
precision ranging to the surface as well as to the top of intermediate cloud layers, which 
is a requirement for the inversion of column CO 2 number density measurements to 
column C02 mixing ratios, has been demonstrated using new hyperfine interpolation 
techniques that takes advantage of the periodicity of the modulation waveforms. This 
approach works well for both BPSK and linear swept-frequency modulation techniques. 
The BPSK technique under investigation has excellent auto-correlation properties while 
possessing a finite bandwidth. A comparison of BPSK and linear swept-frequency is also 
discussed in this paper. These results are extended to include Richardson-Lucy 
deconvolution techniques to extend the resolution of the lidar beyond that implied by 
limit of the bandwidth of the modulation, where it is shown useful for making tree 
canopy measurements. 


1. Introduction 

The ability of lidar to penetrate to the ground through gaps in tree canopies makes it 
especially attractive for generating bare earth terrain models in forested areas. In addition, 
lidar can tell us much about the vegetation itself. Depending on the sensor type and 
collection parameters, a lidar system can either record multiple returns from a single 
transmitted pulse from the top of the tree canopy, intermediate layers of the canopy, as 
well from the ground, or can record all of the laser energy returned in a continuous 
waveform. 

Attempts to use lidar to map forests with profiling laser altimeters first appeared in 
literature in 1980s [1,2] and in the 1990s, a number of efforts validated the use of lidar 
for forestry inventories [3,4,5], Currently, many forestry firms and government agencies 



throughout the world use lidar technology for mapping, inventory, and resource 
management [6, 7, 8]. 

For trees, canopy height is the most common forest characteristic that is of importance. 
Since so much is known about the biophysical attributes of individual tree species, height 
is the driver in estimating wood volume, evaluation of local vegetation productivity, 
estimating biomass, predicting potential yield, commercial value, and fire behavior. Lidar 
has proven to be very successful in estimating tree heights. The measurement of tree 
height for an individual tree is defined as the vertical distance from the base of the tree to 
the topmost terminal bud. The measurement of tree height is estimated from lidar by 
subtracting the highest elevation return from the lowest elevation return closest to the tree. 

For clouds it is extremely important to measure cloud height and thickness especially 
for low level clouds which are a key for cloud climate feedbacks [9, 10]. Measurements 
of the height and thickness of these low level clouds would provide direct information on 
the surface level moisture and boundary layer turbulence since the marine boundary layer 
clouds are generally produced by adiabatic lifting processes of surface moisture. The 
base and top of the clouds represent the lifting condensation level and buoyancy potential 
of surface moist air [11, 12], For high clouds deep convection activities and rainfall 
efficiency, driving factors for high cloud feedbacks, are related to the cloud heights [13]. 
Advanced knowledge of physical properties of boundary layer and high clouds will 
significantly improve model predictions of cloud feedbacks and future climates. 

One problem is that the measurements typically require a vertical-sampling resolution 
of 0.15-0.3 m and a pulse width of 3-10 meters to obtain useful information about trees. 
For clouds, thin cirrus clouds can be as thin as 60 m and boundary layer clouds can be of 
the order of 200-300 m. This limits lidar’s usage in forest and some cloud applications. 
With signal processing, some other lower resolution lidars that were designed for a 
different purpose may be used. In a recent study the Richardson-Lucy algorithm, which is 
normally used to enhance astronomical images [14], was proposed to improve the 
resolution of a high-resolution lidar for the measurement of tree canopies [15]. However, 
this lidar had a sampling resolution of 0.15 m and a pulse width of 2-16 m and the 
improvement was only moderate. For our preliminary focus, we have been working to 
develop a lidar system to measure atmospheric CO 2 for the ASCENDS space mission [16, 
17, 18]. As we have shown in recent flight tests, the vertical-sampling resolution of even 
a low-resolution lidar with a digitized-sampling resolution of 75 m can be enhanced to as 
high as 0.15 m under our processing conditions [19] and the synthetic pulse width can be 
improved from a measured 352 m to less than 5 m [20]. The spatial sample rate 
determines the ability to measure distances from a hard target and the pulse width 
determines the ability to resolve one object from another (canopy vs. ground for instance). 
Without further processing our ability to resolve different objects in the same frame 
would be limited to about 352 m in this case and our ability to measure range would be 
limited about ±37.5 m. So our super-resolution technique is a major advancement for this 
application. What makes this possible is using a near exact interpolation we have 
developed [19] in combination with our new two-stage Richardson-Lucy algorithm [20]. 
This is a two-step approach: enhance vertical resolution first, and, then, deconvolve the 
signals from extended targets such as clouds and forest canopies to obtain super- 
resolution vertical profiles with resolution-enhanced data. Even though the starting 
vertical resolution is crude, the information is embedded in the phase of received signals 



so one can look at many low-resolution profiles of these signals simultaneously to 
reconstruct a high-resolution vertical profile. So far this has been demonstrated on a 
prototype intensity-modulated (IM) continuous wave (CW) lidar designed to measure 
atmospheric CO 2 concentration for the ASCENDS mission, but it could be reformulated 
for low resolution pulse systems such as the Cloud-Aerosol Lidar with Orthogonal 
Polarization (CALIOP) onboard CALIPSO satellite or some other lidar systems. Another 
possible application for the algorithm would be to build a low data rate lidar system 
capable of making atmospheric measurements, where the data rate is low enough to 
downlink the raw data in full and to process the data on the ground, making this a basis 
for a new compressive sensing system. 


2 Super-resolution algorithm description 

2. 1 CW Measurement concept 
The concept of our IM-CW system designed 
for CO 2 integrated path differential absorption 
(IPDA) measurements is illustrated in Fig. 1. 

For the system described in this Fig. 1, 
intensity modulated online and offline seed 
lasers with distinct spectral properties are 
combined using fiber optics and 
simultaneously amplified using a single 
Erbium Doped Fiber Amplifier (EDFA) to 
increase the transmitted power. A small 
fraction of the transmitted beam is picked off 
via an optical tap inside of the EDFA and sent 
to a reference detector for power normalization. 

The backscattered science signals of the online 
and offline wavelengths from the surface as 
well as aerosols and clouds are simultaneously 
collected with a telescope, optically filtered 
with a narrow band optical filter, and detected 
by a single detector. Both the science and 
reference signals are amplified, electronically 
filtered and then digitized for retrievals of 

column CO 2 using IPDA approach. Post processing of the digitized science and 
reference data allows for discrimination between ground and intermediate scatterers using 
the matched filter technique. One also obtains differential absorption power ratios for 
inference of CO 2 column amounts as well as range estimates to the scattering targets. The 
on-line and off-line laser modulation signals can be expressed as 



Fig. 1. Baseline instrument block diagram. 
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where m is the modulation depth with a value between 0 and 1. The quantity §0) is the 
repeating modulation waveform where -l<£(r)<l and (£(r)} = 0. In general, the on- and 

off-line waveforms are chosen to be orthogonal to each other to avoid interference 
between channels. 


Given the modulation waveforms presented in Eq.l, the instantaneous received lidar 
off-line and on-line powers (p'h)and p’ ( t ) , respectively) from a single hard scatterer at 

range r are given by 
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where t is the total one-way column optical depth for the off-line measurement that is 
determined by the gas absorption from molecules other than CO 2 , and t ' is the one-way 
column optical depth resulting from CO 2 absorption only, [3 is the backscatter coefficient 
(km 'sr 1 ), S is the extinction to backscatter ratio (sr 1 ), Kk is a function of the kth 
scatterer’ s reflectivity, and P' off and P' m are the average transmitted off-line and on-line 

powers over one period of the modulation waveform, respectively. 

When multiple targets including the surface exist in the lidar path length, each target 
would generate similar received lidar signals like those in Eq. 2. The signals from these 
multiple targets are combined at the detector and converted to an electronic signal s(t). 
For an AC coupled receiver subsystem, we have, 

s (t) = X[ c u<„ (* - 2r k I c) + C 2k m$ off ( t - 2 r k / c)] (3) 

k 

where we sum over all k scatterers and where, 
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where K’ is a constant. These Cik and C 2 k returns can be uniquely discriminated from 
other returns from different scatterers using the matched filter with the transmitted 

waveform. These are also proportional to the average received optical power ( P* ff k and 
p l , k ) from the kth scatterer. For a ground target, solving for t g ' gives 
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where C lg and C 2g are the ground returns. Once C lg and C 2g are determined, the column 

optical depth for CO 2 can be found using Eq. 5. In general, this is done by cross 
correlating the reference waveform with the return signal s(t), i.e., by the matched filter 
technique as mentioned before. 

2.2 Multi-swept Frequency modulation 



The main advantage of the multi-sweep approach is one may achieve a high level of 
orthogonality between channels while using simultaneous channels that are very close in 
frequency. The advantage is that one may make simultaneous on-line/off-line 

measurements with virtually no crosstalk. 

The modulation waveform for the kth orthogonal channel takes the form 

&(0 = cos (&( 0 )> ( 6 ) 

where 

<p k {t) = 2n\‘j k {t')dt' , (7) 

where fk is the frequency for the k t h channel defined by 

f k (t)=f 0k +fs„ee P { t -T M'/r)), (8) 

where T is the sweep period, f 0k is the start frequency, f meep is the sweep function, and 
int(?/7’) is the integer part of t/T . Let 

^ee P (t)=2Kf o f swe Jt')dt'. (9) 


We now construct a continuous phase function for each channel such that 

^(0 = 2^/ ot t+^(r)int(t/7’)+^(f-r int(f/r)). (10) 

We digitize these waveforms at a sample rat ef s , a sample time of A? = 1 / f s and with a 
frame size of N points that includes M sweeps such that the total time for M sweeps is 
MT = NAt . We define a circular correlation such that 
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where ref is the reference waveform, data is the data collected either from the reference 
or science detector, and DFT is the digital Fourier transform. For our data processing, ref 
will take the form of exp(/ <(>{) because our correlation is performed in quadrature. The 

channel frequencies f 0k are chosen to make each channel orthogonal with every other 
channel such that 


/?(exp(/ (j)j ), cos(d t )j = 0 , j^k 


( 12 ) 



where / = V-T . Since each channel is orthogonal, the constants C\ and C 2 will then be 
proportional to the peak of the cross correlation of the return signal with the reference 
kernel for channels 1 and 2 respectively and the range is determined by the delay in the 
correlation where distance is r = c 8tl 2 where c is the speed of light and* is the time 
delay. 

The general swept case is closely related to the unswept case [21] in that the choice of 
frequencies is similar if we use the average frequency over a single sweep. That is to say 
one may use average frequency of 
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with the results of the unswept case as a starting point to guess what the orthogonal start 
frequencies are. Though related, the swept case is somewhat more complicated in the 
choice of allowable start frequencies. The role of necessity in choosing start frequencies 
can be presented as 
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where K is the number of channels and nj is an integer. The autocorrelation function for 
each channel will have exactly M peaks. Generally speaking, if the number of channels is 
greater than 2, the number of sweeps must be greater than the number of channels (M>K). 
For the two-channel case we must have at least 2 sweeps. The integer^ is chosen to make 
the autocorrelation for the first channel have the same peak height for each sweep. Once 
this is done, n j for the other channels are chosen to make each channel orthogonal with 

every other and the autocorrelation function for each channel has equal peak heights. In a 
four channel setup, for instance, this can usually be done quickly by evaluating the cross 
correlation between channels and autocorrelation functions. 

For the linear sweep case we choose 

L,W = 4fpO<Kr, (15) 

where A f is the sweep bandwidth. This is a well-known case from radar and will result in 
a range resolution of 
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where c is the speed of light. The unambiguous range is given by 
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The frequency as a function of time for each channel is 


where int [t IT) is the integer part of t IT . The phase is 


In the continuum limit it can be shown the autocorrelation function shown in Fig. 2 is 
given by 


In the above integration we ignore the 
lower sideband because that is mostly filtered °- 5 
out in the integration. 0 .4 

As an example we take the special case . 
where we have a sample rate of 2000 kHz, a f 

sweep bandwidth of 500 kHz, a frame size of < 0 2 

4096, and a sweep period of 512/2000 ms „.i 
which results in M=8 sweeps per frame. We 


Care must be taken to compute these as accurately as possible or choose parameters that 
result in start frequencies that are near integers to minimize round off error. Cross 
correlation crosstalk is 0 to within round off error in that case and between all other 
channels as well. 
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chose 4 channels and the parameters chosen 
that result in orthogonal sweeps are «i=1450, 
«2=18, a 3 =30, and m=52. This results in start 
frequencies of approximately 104.004 kHz, 
108.398 kHz, 111.328 kHz, and 116.699 kHz. 
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Fig. 2. Autocorrelation function using linear 
swept frequency modulation. 


2.3 Binary Phase Shift Keying (BPSK) 

The BPSK approach discussed here is where an IM sine wave carrier is modulated by a 
pseudo-noise (PN) code in the form of a maximum length (ML) sequence as in Eq. 21. 
To effectively represent each bit of ML-sequence in the transmitted waveforms, and 
allow flexibility in the unambiguous range for a particular code length and sample rate, 
we oversample the ML-sequence by an integer M such that each code bit of the ML 
sequence is represented by M points. Let z„ = z(n) be the original ML-sequence and 

z„=z(n)be the oversampled ML-sequence such that z(n) = z(int((n-l)/M)+l) , where 

int(x) is a function that represents the integer part of x. The on and offline modulation 
waveform is given by, 


&(0 = &(«■ 1 f s ) = {2Z(n)-\)co%{2nn /,//;), (21) 

where f on is the online carrier frequency, f off is the offline carrier frequency, and f s is the 

sample rate. We could have more channels for different applications. A previous 
publication [21] describes a procedure for choosing orthogonal swept frequency signals 
where the starting frequency of the sweep is the desired parameter. We modify that 
slightly by choosing the center frequency as the desired parameter for a BPSK signal. Let 
us suppose we have exactly P repeats of the ML sequence in a frame of N points where 
the period of each repeat is T. We claim the center frequency of K orthogonal channels 
will take the form 
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Not every choice of nk will lead to orthogonal frequencies but the ones that are 
orthogonal will take that form. This form works best in situations where there are no 
multiple scatterers or where only ranging is needed as under clear skies. But a small flat 
background offset proportional to the inverse of the PN code length could exist. This 
does not present a problem unless there are multiple scatterers, where it could cause a 
small bias error in the peak height due to interference with other scatters. Choosing a very 
long PN code can minimize this effect. This reference waveform is given by 

r k(n) = {2Z(n)-l)exp(27tinf on / f s ). (23) 


This improved form can be used to both clear and cloudy sky cases and has the advantage 
of a perfectly 0 background offset so there is no interference between scatters and no bias 
in the peak heights. In this case the reference signal is 


r '* («) = Z(n)e\p(2jcinf k / f s ) . 


(24) 
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Fig. 3. Autocorrelation function using the BPSK 
modulation is completely lacking in sidelobes or 
other artifacts, which is the key advantage for 
using BPSK modulation. 


Demodulation is accomplished simply by 
taking the magnitude of the matched fdter 
correlation described by Eq. 1 1 . 

As an example of a modulation 
technique that is suitable for a band pass 
fdtered LAS system presented in Fig. 1, 
we take the special case where P=16 (16 
ML sequence repeats) and M=4 (4 samples 
per code bit) and choose a 2 MHz sample 
rate. We generate a 7 th order (or 127 
length) ML-sequence using the recurrence 
relation z n+1 = z„ ®z n+6 with the seed {1,0, 1, 
0, 1, 1, 1}, where© is the “exclusive or”. 


An example of a 6-channel system that is 
mutually orthogonal using a frame size of 8128 
(16, 127 point ML sequence repeats that were 
oversampled by a factor of 4) would be to 
choose center frequencies of {225125/508 kHz, 
56875/127 kHz, 230125/508 kHz, 232625/508 
kHz, 117625/254 kHz, 237625/508 kHz}. 
These are displayed as integer ratios so they 
may be computed in a more exact fashion to 
produce perfectly orthogonal BPSK waveforms. 
A plot of the autocorrelation function (derived 
from Eq. 11 and Eq. 24) is shown in Fig. 3, 
which is ideal in that it has no sidelobes or 
other artifact and is the synthetic pulse derived 
from the BPSK modulation. 
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Fig. 4. Frequency spectrum of correlation 
before reordering shows spectrum is 
formed from a comb as we see from the 
blowup of the of the spectrum which 
shows a 25 point spacing between the teeth 
of the comb. 


2.4 Fourier Transform Reordering (FTR) interpolation of periodic waveforms 
Given that the DFT of the correlation forms a comb in the case of a repeating waveform 
and that the spacing between the teeth of the comb is equal to the number of repeats, we 
have shown that if one reorders the DFT array elements such that the spacing between the 
teeth is 0, one can effectively transform a periodic pulse train into a highly interpolated 
single pulse [19]. However, this must be done about a point equal to the center frequency 
of the spectrum plus the Nyquist rate such that all the teeth to the left of this point are 
moved to the left and all the teeth to the right of this point are moved to the right. This 
works because there is no mirror image in the frequency domain and it will therefore only 
have one unique spectrum since we are using the complex quadrature as the reference. 
The only requirement for this to work is that the signal must be sampled higher than the 
Nyquist rate. The reordered array is 

f\n) = fln 0 +N s (n-l)], 

1 <n<N /2 + N -1, 

p c 7 

f\n + N-N p /2) = f[N s (n-l + N p /2) + n 0 ], (25) 

N <n< N 12, 

c p 

f\n) = 0, otherwise. 
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Fig. 5. Frequency spectrum after reordering with 
the comb tooth spacing reordered to 0 as is shown 
in the blowup. 

Here /’ is the new spectrum, / is the old 
spectrum, no is the position of the first tooth 
of the comb, N s is the number of 
sweeps/repeats, N is the total number of 
points, 7Vp is the number of points in one 
sweep/repeat, and N c is the center frequency 
of the waveform in bins. 

Alternatively we can do a cyclic rotation in the frequency domain such that the center 
frequency is mapped to 0. Then we can reorder the array such that all the teeth left of 
center are moved to the left and all teeth to the right of center are moved to the right. In 


Fig. 6. Result of correlation before 
reordering shows periodic pulse profile with 
crude sampling. 


this case one would use Eq. 25 with N e =0. Both of the aforementioned methods work - 
the only difference is the phase term in the resulting correlation, which is removed when 
the magnitude is taken. The level of interpolation will be equal to the number of repeats 
in the DFT. 

As an example we take a linear sweep case with 10000 points total with 25 sweeps. Fig. 4 
shows a plot of the amplitude of the DFT of the correlation as described in Eq. 1 1 and 
shows how using the complex quadrature gives only one unique spectrum without a 
mirror image, thereby extending the Nyquist rate by a factor of 2. This is one of the many 
advantages to performing the correlation in quadrature. The blowup in Fig. 4 shows this 
spectrum is actually a comb. The tooth spacing in this case is 25 because we have 25 
sweeps. Fig. 5 shows the spectrum after reordering with a comb tooth spacing of 0, which 
shows a more continuous curve. 

Fig. 6a shows the resulting correlation before reordering. Since we are correlating 
multiple sweeps at once the result is multiple pulses. One can see from Fig. 6b the detail 
of each pulse shows crude sampling. Fig. 7a shows the resulting correlation after 

reordering results in a single pulse. The 
detail of this pulse shown in Fig. 7b 
exhibits fine sampling resulting from a 
factor of 25 improvement in the apparent 
sample rate. 

The theoretical explanation for this is the 
DFT of a function that is periodic within 
the range of the DFT forms a Kronecker 
comb. By changing the comb spacing to 0, 
the resulting inverse Fourier transform 
becomes non-periodic within the DFT 
range. Since the values outside of the comb 
are 0, this operation performs a type of 
Fourier interpolation. However, unlike 
normal Fourier interpolation we don’t 
introduce extra points or 0s in the spectral 
domain - we just reorder the ones that 
already exist, which makes this technique 
very efficient. This technique [19] also 
saves tremendous computational time in 
order to increase ranging precision and 
accuracy compared to the time domain 
oversampled by the number of repeats, which curve fitting and interpolation technique 
in this case is 25. used previously [17] [18]. 

2.5. Super-resolution by FTR interpolation and Richardson-Fucy (RF) deconvolution 
After we take the magnitude of the resulting interpolated correlation obtained from Eqs. 
1 1 and 25 we may further enhance this by implementing an accelerated RF algorithm 
[22] to deconvolute lidar returns from extended backscatteres for tree height and 
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Fig. 7. Result of correlation after reordering 
shows crudely sampled periodic pulse train is 
converted into a single pulse that is 


thickness estimations. The basic assumption used by RL is the measured image is given 
by 


A-y®p + ri 


(26) 


where A is the measured image, y is the Point Spread Function (PSF), p is the 
undistorted image, r\ is the noise, and ® is the convolution operator. The RL algorithm 
estimates the undistorted image iteratively by calculating N iterations of 
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where p„ is the nth estimate of the undistorted 
image and ® is the correlation operator. Both 
the PSF and image are assumed to be positive 
with RL. Because we are using such a high 
degree of interpolation, we modify the 
algorithm further to second order by iterating 
with a new PSF derived from the original 
PSF such that we also calculate N iterations 
of 
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where y/\=y/ . Then we refine the result 
further by calculating W iterations of 


Fig. 8. Results of FTR interpolation and the 
two-stage RL deconvolution on the ground 
calibration measurement (a) and the zoomed in 
version (b). Using 100 first order iterations and 
30 second order iterations the sample resolution 
goes from 75 m to .1875 m and the measured 
half height pulse width changes from 352 m to 
4.40 m - almost 2 orders of magnitude. Data 
has been centered about 15 km. 
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In this way it is possible to accelerate 
convergence by an order of magnitude or 
more beyond what is described by Biggs et. 
al. [22], 

This algorithm was recently tested during a 
series of airborne test flights for a newly developed IM-CW lidar system known as the 
ASCENDS CarbonHawk Experiment Simulator (ACES) [23]. The demonstration of 
these concepts was carried out where the instrument was operated at ~ 10 W avg. total 
power for three orthogonal transmitted channels and high altitude (~ 9-11.5 km). These 
three orthogonal waveforms were transmitted at three different wavelengths (1 on-line 
and 2 off-line wavelengths in the 1.571 pm CO 2 absorption band) where each transmitted 
wavelength had approximately 3.8 W, 3.5 W, and 2.7 W respectively of average power. 
Using the BPSK modulation described in section 2.2, the instrument operated at a sample 
rate of 2 MHz and a modulation code bitrate of 500 kbps, with 400 x 508 point repeats 



making a total frame length of 203200 [24], which gives us a sampling resolution of 75 
meters and a theoretical pulse resolution of 300 m. Using FTR interpolation we transform 
a periodic 400-peak pulse train to a single pulse interpolated/sampled at 75/400 = 0.1875 
m. 

For the first test we did a ground 
experiment, where the laser was fed back to 
the receiver at close distance at low power. 

By taking the measured ground pulse and 
applying 100 1st order iterations using Eq. 

27 and Eq. 28, then another 30 2nd order 
iterations using Eq. 29, we transform the 
pulse from a measured half height width of 
352 m to 4.40 m, making the width 
resolution almost 2 orders of magnitude 
better as in Fig. 8. 

The PSF was obtained by taking the 
theoretical ML-sequence, autocorrelation 
function, then filtering it with a periodic 
Gaussian filter [25]. Through a 
combination of instrument characterization 
and filtering to obtain the minimum pulse 
width from a hard ground return while still 
being able to see cloud and tree returns, we 
were able to obtain consistent results. 

We conducted many flight tests over 
ocean, land, and clouds. Fig. 9 shows a 
case where we transmitted through thin 
cirrus clouds and reflected from a hard 
surface. Despite the crude sampling 
resolution we are able to discern features 
well below the original 75 meter sampling 
and 352 m pulse width. One interesting 
feature is the change in relative pulse 
height. This can be explained by the 
algorithms tendency to preserve the area 
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Fig. 10. RL enhanced pulse from tall trees 
shows tree height of about 30 m. 


Fig. 9. Transmission through thin cirrus clouds 
(a) shows FTR interpolation combined with RL 
deconvolution is able to detect the thickness of 
thin cirrus clouds (b), yet the half height pulse 
width from the ground is also measured at about 
4.91 m (c), which is only slightly wider than the 
4.40 m enhanced pulse measured in ground 
testing. 

under the pulse. Since the ground pulse 
received from hard surface is compared to the 
thickness of the clouds, as one enhances the 
resolution the ground pulse becomes 
narrower and taller in comparison to the 
cloud pulse. 

As a comparison we enhanced a ground 
return from forest trees, which is shown in 


Fig. 10. Here we see an estimated tree height of about 30 m, which is consistent with the 
type of trees found in southern Virginia. 

2.6 Extentions to pulse lidar 

Thus far we have discussed CW lidar as has been implemented in the Langley 
ASCENDS effort. However, there is no reason this theory couldn’t be modified to cover 
pulse lidar systems such as those used on CALIPSO and other programs, and since this 
type of lidars is more common there is a huge amount of data one could use to extract 
science from. Technically the difference is fairly simple to implement. Instead of using a 
reordered FFT of the correlation to do the interpolation, one would use a reordered FFT 
of a periodic pulse train multiplied by the appropriate comb filter. The Richardson-Lucy 
part of it would be identical. If awarded, we would extend this theory and our algorithm 
used in ASCENDS to include the most common pulse lidar systems and process available 
data such as CALIPSO data, Langley High Spectral Resolution Lidar (HSRL)[26] pulse 
data, etc. 

3. Conclusion 

These results show it is possible to increase both the sample rate and fundamental 
resolution from a CW ranging lidar where the signal is sampled greater than the Nyquist 
rate. Using the Fourier transform reordering interpolation technique to increase the 
sampling resolution and by iteratively applying the modified Richardson-Lucy technique 
presented here, it is possible to increase the resolution by almost two orders of magnitude 
in as has been demonstrated here. This makes it possible to measure tree canopy heights, 
and also opens the possibility of a compressive sensing lidar. 

References 

1. R. Nelson, W. Krabill, & G. MacLean, Determining forest canopy characteristics 
using airborne laser data. Remote Sensing of Environment, 15, 201-212 (1984). 

2. R. Nelson, W. Krabill, & J. Tonelli, Estimating forest biomass and volume using 
airborne laser data. Remote sensing of environment, 24 , 247-267 (1988). 

3. E. Naesset, Estimating timber volume of forest stands using airborne laser scanner 
data, Remote Sensing of Environment, 61, 246-253(1997). 

4. M. A. Lefsky, W. B. Cohen, S. A. Acker, G. G. Parker, T. A. Spies, D. Harding, 
Lidar remote sensing of the canopy structure and biophysical properties of Douglas- 
fir western hemlock forests, Remote sensing of environment, 70, 339-361(1999). 

5. J. B. Blair, D. L. Rabine, M. A. Hofton, The Laser Vegetation Imaging Sensor: a 
medium-altitude, digitisation-only, airborne laser altimeter for mapping vegetation 
and topography, ISPRS Journal of Photogrammetry and Remote Sensing, 54 , 115- 
122(1999). 

6. M. A. Wulder and Seemann, D. Forest inventory height update through the 
integration of LIDAR data with segmented landsat imagery. Can. J. Remote Sens, 
29, 536-543(2003) 



7. E. Naesset, T. Gobakken, J. Holmgren, H. Hyyppa, J. Hyyppa, M. Maltamo, U. 
Soderman, Laser scanning of forest resources: the Nordic experience, Scandinavian 
Journal of Forest Research, 19, 482-499(2004). 

8. H. E. Andersen, R. J. McGaughey, S. E. Reutebuch, Estimating forest canopy fuel 
parameters using LIDAR data. Remote sensing of Environment, 94, 441-449(2005). 

9 Steven C. Sherwood, Sandrine Bony, and Jean-Louis Dufresne, Spread in model 
climate sensitivity traced to atmospheric convective mixing, Nature, 55, 37-42, 
2014. 

10 Bony, S. & Dufresne, J. L. Marine boundary layer clouds at the heart of tropical 
cloud feedback uncertainties in climate models. Geophys. Res. Lett. 32, L20806, 
2005. 

11 Lin, B., P. Minnis, and A. Fan, Cloud liquid water amount variations with 
temperature observed during SHEBA experiment, J. Geophys. Res., 108 (D14), 
4427, doi:l 0.1 029/2002 JD002851, 2003. 

12 Lin, B., P. Minnis, T.F. Fan, Y. Hu, and W. Sun, Radiation characteristics of low 

and high clouds in different oceanic regions observed by CERES and MODIS, the 
International Journal of Remote Sensing, 31, 6473-6492, DOI: 

10.1080/01431160903548005, 2010. 

13 Lin, B., B.A. Wielicki, P. Minnis, L. Chambers, K.-M. Xu, Y. Hu, and A. Fan, The 
effect of environmental conditions on tropical deep convective systems observed 
from the TRMM satellite, J. Climate, 19, 5745-5761, 2006. 

14. L. B. Lucy, An iterative technique for the rectication of observed images, The 
Astronomical Journal, 79, 745-754(1974). 

15. J. Wu, J. A. N. van Aardt, G. P. Asner, A comparison of signal deconvolution 
algorithms based on small-footprint LiDAR waveform simulation, Geoscience and 
Remote Sensing, IEEE Transactions on, 49, 2402-2414(2011). 

16. NRC, Earth Science and Applications from Space: National Imperatives for the 
Next Decade and Beyond, The National Academies Press, Washington, D.C., 2007 

17. B. Lin, S. Ismail, F. W. Harrison, E. V. Browell, A. R. Nehrir, J. Dobler, B. Moore, 
T. Refaat, S. A. Kooi, “Modeling of intensity-modulated continuous- wave laser 
absorption spectrometer systems for atmospheric CO 2 column measurements,” Appl. 
Opt., 50(29), 7062-7077 (2013). 

18. J. T. Dobler, F. W. Harrison, E. V. Browell, B. Lin, D. McGregor, S. Kooi, Y. Choi, 
and S. Ismail, “Atmospheric CO 2 column measurements with an airborne intensity- 
modulated continuous wave 1.57 pm fiber laser lidar,” Appl. Opt., 52(12), 2874- 
2892 (2013). 

19. J. F. Campbell, B. Lin, A. R. Nehrir, F. W. Harrison, M. D. Obland, High 
Resolution CW Lidar Altimetry using Repeating Intensity Modulated Waveforms 
and Fourier Transform Reordering, Opt. Lett., 39, 6078 (2014). 

20. J. F. Campbell, B. Lin, A. R. Nehrir, F. W. Harrison, M. D. Obland, Super- 
resolution technique for CW lidar using Fourier transform reordering and 
Richardson-Lucy deconvolution, 24, 6981 (2014) 

21. J. F. Campbell, Nonlinear swept frequency technique for CO 2 measurements using 
a CW laser system, Appl. Opt., 52, 3100 (2013). 



22. D. S. C. Biggs and M. Andrews, Acceleration of Iterative Image Restoration 
Algorithms, Appl. Opt., 36, 1766(1997). 

23. Obland, M., C. Antill, E. V. Browell, J. Campbell, S. Chen, C. Cleckner, M. 
Dijoseph, F. Harrison, S. Ismail, B. Lin, B. Meadows, C. Mills, A. Nehrir, A. Notari, 
N. Prasad, S. Kooi, N. Vitullo, J. Dobler, J. Bender, N. Blume, M. Braun, S. Homey, 
D. McGregor, M. Neal, M. Shure, T. Zaccheo, B. Moore, S. Crowell, P. Rayner, W. 
Welch, 2013 American Geophysical Union Fall Meeting, San Francisco, CA, 
December 9-13, 2013. 

24. J. F. Campbell, B. Lin, A. R. Nehrir, F. W. Harrison, M. D. Obland, Binary Phase 
Shift Keying on Orthogonal Carriers for Multi-Channel CO 2 Absorption 
Measurements in the Presence of Thin Clouds, Opt. Exp., 22 , A1634 (2014). 

25. J. F. Campbell, B. Lin, A. R. Nehrir, Advanced sine wave modulation of continuous 
wave laser system for atmospheric CO 2 differential absorption measurements, Appl. 
Opt., 53, pp. 816-829(2014). 

26. J. W. Hair, C. A. Hostetler, A. L. Cook, D. B. Harper, R. A. Ferrare, T.L. Mack, W. 
Welch, L. R. Izquierdo, and F. E. Hovis, Airborne High Spectral Resolution Lidar 
for profiling aerosol optical properties, Appl. Opt., 47, pp. 6734-6752 



